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We investigate the 0(a s v 2 ) correction to the e + e~ — > J/ip + n c process in the nonrelativistic QCD 
(NRQCD) factorization approach. Within some reasonable choices of the relative order-v 2 NRQCD 
matrix elements, we find that including this new ingredient of correction only mildly enhances the 
existing NRQCD predictions. We have also deduced the asymptotic expressions for the 0(a s v 2 ) 
short-distance coefficients, and reconfirm the early speculation that at next-to-leading order in a s , 
the double logarithm of type In 2 (s/m 2 ) appearing in various NRQCD short-distance coefficients is 
always associated with the helicity-suppressed channels. 
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I. INTRODUCTION 

One of the most intensively studied hard exclusive re- 
actions in recent years is perhaps the double-charmonium 
production process e + e~ — » J/iJj + i] c at the B fac- 
tory energy y/s — 10.58 GeV [1]. It was initially mea- 
sured by the Belle experiment in 2002 with a[e + e~ —> 
J/tp+r) c ]xB>A = 33±g±9 fb [2], where B> A is the branch- 
ing ratio of r\ c into 4 or more charged tracks. The first 
theoretical predictions [3-5], built on the lowest-order 
(LO) calculation in the nonrelativistic QCD (NRQCD) 
factorization approach [6], were scattered in the range 
2.3 — 5.5 fb, almost one order of magnitude smaller than 
the Belle data. Later Belle Collaboration refined their 
measurement and gave a[J/ip+rj c ] x2? >2 = 25.6±2.8±3.4 
fb [7], where S>2 denotes the branching fraction for the 
rjc into more than two charged tracks. In 2005, BaBar 
Collaboration also measured the same observable and ob- 
tained 17.6±2.8i2:? fb [8]. 

The disquieting discrepancy between experiment and 
the LO NRQCD predictions has spurred a great amount 
of theoretical endeavors in the following years. Roughly 
speaking, most works can be divided into two major cate- 
gories, either based on the light-cone factorization [9-12], 
or based on the NRQCD factorization approach [13-16] 
(for the investigations from other theoretical approaches, 
see Refs. [17-20]). 

One crucial step toward alleviating the tension between 
data and theory is the discovery of the positive and sig- 
nificant 0(a s ) correction to the e + e~ — > J/ip + i] c pro- 
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cess [13, 14]. This ncxt-to- leading order (NLO) perturba- 
tive calculation was performed in the NRQCD factoriza- 
tion framework (a similarly large positive 0(a s ) correc- 
tion has also later been found for the e + e~ — > J/ip + Xco 
process [21, 22]). In contrast, due to some long-standing 
theoretical difficulty inherent to the helicity-suppressed 
process, by far no one has successfully conducted the cor- 
responding 0(a s ) correction to this process in the light- 
cone approach. In a sense, for the e + e~ — > J/ip + rj c 
process, NRQCD approach seems more systematic and 
maneuverable than the light-cone approach. 

The relative 0(v 2 ) correction to e + e~ — > J/ip + rj c has 
also been addressed [3, 15, 16], where v denotes the typ- 
ical velocity of the c quark in a charmonium. Notwith- 
standing the large uncertainty inherent to relativistic cor- 
rection, it was believed that [15, 16], including both NLO 
perturbative and (a partial resummation of ) relativistic 
corrections, one may achieve the reasonable agreement, 
albeit with large uncertainties, between the NRQCD pre- 
diction and B factory data. 

The goal of this work is to address the 0(a s v 2 ) cor- 
rection to e + e~ — > J/t/j + rjc in the NRQCD factorization 
approach. It is curious to examine its phenomenological 
impact. On the other hand, thus far there are only very 
few basic quarkonium decay processes whose 0(a s v 2 ) 
corrections have been calculated, e.g. J/tp —> e + e~ [23], 
and rj c — > 77 [24, 25]. Therefore, it is also theoretically 
interesting to know the 0(a s v 2 ) effect for the exclusive 
quarkonium production process in the first time. 

The rest of the paper is structured as follows. In 
Sees. II and III, we express the product rate for e + e~ — > 
J/ip + rjc in terms of the J/ip + rj c electromagnetic 
(EM) form factor, and present the NRQCD factorization 
formulas for both quantities, accurate through relative 
order-ii 2 . In Sec. IV, we list the tree-level short-distance 
coefficients through relative 0(v 2 ). In Sec. V, we first 



2 



sketch some key technical steps about the NLO perturba- 
tive calculations, then present the asymptotic expressions 
for the matching coefficients. We devote Sec. VI to ex- 
ploring the phenomenological impact of our new 0(a s v 2 ) 
correction on the B factory measurement. Finally, we 
summarize in Sec. VII. 



II. J/V + ri c EM FORM FACTOR 

Suppose we work in the e~ and e + centcr-of-mass 
frame with invariant mass of ^/s. Let P\ (A) denote the 
momentum (helicity) of the J/tp, and P 2 the momentum 
of the r\ Cl respectively. This process simply probes the 
J/ip + r/c electromagnetic form factor in the timelike re- 
gion, referred to as G(s) hereafter: 

(J/<P(P 1 ,\)+r ]c (P 2 )\j£ M \0) = iG{s)e^P lv P 2p el{X)ll) 

where J^ M is the electromagnetic current. The tensor 
structure specified in (1) is uniquely dictated by the 
Lorentz and parity invariance. As a result, the outgo- 
ing J/ip must be transversely polarized, i.e., A = ±1. 
The cross section can be expressed as 



a\e + e 



J/ip + Vc] 



\G{s)\\ (2) 



where |P| signifies the magnitude of the 3-momentum 
carried by the J/ip (rj c ) in the center-of-mass frame. The 
cubic power of |P| is reminiscent of the fact that J/tp and 
rjc are in the relative P-wave orbital state. 

It is worth recalling the asymptotic behavior of (2) 
for this helicity-flipped process. The helicity selection 
rule [26] dictates that G(s) <~ 1/s 2 as y/s >• m c [10], 
hence a[J/ip + r] c ] ~ 1/s 4 . Evidently, the charm quark 
mass may serve as the agent of violating the hadron he- 
licity conservation. As will be examined in Sec. V, this 
power-law scaling is subject to double-logarithmic mod- 
ification once beyond LO in a s . 



III. NRQCD FACTORIZATION FORMULA FOR 
J/ip + »7c PRODUCTION RATE 

According to NRQCD factorization formula, the J/tp+ 
7y c EM form factor in (1) can be factorized as 



G(s) = ^Mj^M^J/^a ■ e X \0)(Vc\^x\0) 

[cq + C2,i(v 2 ) J/4 , + C2^(v 2 ) Vc H ] , (3) 

where Co and a are the corresponding short-distance co- 
efficients. We have adopted relativistic normalization for 
the quarkonia states appearing in the left side, while 
using the nonrelativistic normalization for those in the 
NRQCD matrix elements. For simplicity, we have in- 
troduced the following dimensionless ratios of NRQCD 



matrix elements to characterize the 0{v 2 ) corrections: 
, 2 . (J/jj(XMH-0) 2 ^-e(X) X \O) 



ml (J/VWto- . 6 (A) X |0) 



_ ( Vc \^(-ir>r x \o) 

\ v )vc - 



ml (»7 C 1^x10) 



(4b) 



where ip^I^x = V^Dx — (D^)^. 

Substituting Eq. (3) into (2), one can decompose the 
cross section into the O(v ) and 0(v 2 ) pieces: 

a[e + e~ -> J/^Jj + r/ c ] = a Q + a 2 + O(a v 4 ), (5) 



where 



87raV(l-4r) 3 / 2 2 
= 5 (Ui)j/4>{Ui)riJco\ , (6a) 



47ra 2 m 2 (l-4r) 3 /2 



(Oi) W <Oi>, c (6b) 



{(T^" |co|2 + 4Re[coC2 ' l] ) (w2>J/,A 



1 - 10r 



+ V7^ |co1 +4Re[c ° C2 < 2] ] ( ^ r 



In deriving (6), we have employed the Gremm-Kapustin 
relation [27] M% w 4m 2 (l + (v 2 ) H ) to eliminate the 
explicit occurrences of Mj/^ and M Va . For notational 
brevity, we have introduced the following symbols: r = 

ip x ) JH , = \m\^cr ■ CX |0)| 2 , and (O^ = 

|fa# + x|o>| 2 . 

It is convenient to organize the short-distance coeffi- 
cients Ci in power series of the strong coupling constant, 
i.e., Ci = c- ' + — c- 1 ^ + • • • . Accordingly, we may decom- 



pose the cross section Ui into + a\ L> (i = 0, 2). Our 
primary goal in this work is to calculate . 



IV. TREE-LEVEL NRQCD SHORT-DISTANCE 
COEFFICIENTS 



It is straightforward to employ the perturbative match- 
ing method to determine the short-distance coefficients, 
by replacing the J/V> and r\ c states with the free cc( 3 s[^) 
and cc( 1 5^ 1) ) pairs, and enforcing that both perturbative 
QCD and NRQCD calculations in (3) yield the same an- 
swer, order by order in a 8 . 

There are in total 4 diagrams at LO in a s , one of which 
is depicted in Fig. 1(a). The tree- level short-distance co- 
efficients through 0(v 2 ) have been available long ago [3]. 
Here we list their values at D = 4 — 2e spacetime dimen- 
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then apply the covariant spin projector [31] to enforce 
two cc pairs to form the spin-triplet /singlet, color-singlet 
states, with the Dirac and color traces handled by Feyn- 
Calc [32]. 

We then expand the amplitude in powers of the quark 
relative momenta, qi, up to the quadratic order. We then 
make the following substitution to project out the iS-wave 
states: 



a v 

QiQi 



D — 



T (-<r 



p2 , 



(8) 



FIG. 1: One sample LO diagram and five sample NLO dia- 
grams that contribute to 7* — > J/ip + r/ c . 
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(7a) 
, (7b) 
,(7c) 



where e c = § is the electric charge of the charm quark, 



N c — 3 is the number of colors and Cf — 



N- 



2N C 



V. NLO PERTURB ATIVE NRQCD 
SHORT-DISTANCE COEFFICIENTS 

In this section, we first sketch some technical issues 
about the NLO perturbative calculations, then present 
the analytic expressions of the NLO corrections to the 
matching coefficients Co and C2. Technical details will be 
expounded in a long write-up [28] . 



A. Description of the calculation 

The main technical difficulty is to identify the relativis- 
ts effects in the quark amplitude 7* — > cc(P\, 3 Si) + 

cc(P2, 1 Sq 1 ^) through NLO in a s . The momenta of c and 
c inside each pair are assigned as pi = ^Pi + qi, pi — 
\Pi—qi for % = 1,2. The total momentum and relative 
momentum qi are chosen to be orthogonal. In the rest 
frame of each cc pair, the total and relative 4-momenta 
read Pf 1 — (2E qi ,0), gf = (0, qi), respectively, where 

We first employ the Mathematica package Fey- 
nArts [29, 30] to generate Feynman diagrams and am- 
plitudes for the process 7* — > c{pi)c{pi) + c(p2)c(p2) to 
NLO in a s . In total there are 20 two-point, 20 three- 
point, 18 four-point, and 6 five-point one-loop diagrams, 
some of which have been illustrated in Fig. !(&)-(/). We 



for i — 1,2, and q 2 is understood to be defined in the 
rest frame of each cc pair. 

In conventional matching procedure, one expands the 
relative momentum qi only after completing the loop in- 
tegrations in the QCD amplitude, which is a daunting 
task in our case since the entanglement of three dis- 
parate scales, y/s, m c , q^, in a loop integral. In this 
work, we employ a much simpler shortcut suggested by 
the method of region [33], i.e., making expansion in qi 
prior to carrying out the loop integration. This amounts 
to directly extracting the NRQCD short-distance coeffi- 
cients, i.e., the contributions solely arising from the hard 
region (A: 2 > m 2 c ). Consequently, we will no longer be dis- 
tracted by the effects from the low-energy regions such 
as the potential (k° ~ mv 2 , |k| ~ mv) region. 

To proceed, we use the Mathematica packages 
Fire [34] and Apart [35] to reduce the general higher- 
point one-loop tensor integrals into a set of master in- 
tegrals. As a bonus of having expanded the integrand 
in powers of qi, and utilized a trick of rescaling the pair 
momentum Pi to make some hidden relativistic effects ex- 
plicit, it turns out that all the required master integrals 
are just the usual 1-, 2- and 3-point scalar integrals, all 
of which can be found in the Appendix of Ref. [14]. 

We adopt the dimensional regularization to regularize 
both UV and IR singularities, with spacetime dimension 
D — 4 — 2e. We use 't Hooft-Veltman scheme to handle 
75 [36, 37]. After summing the contributions from all the 
diagrams, and incorporating mass and coupling constant 
renormalization, we find that the ultimate NLO expres- 
sions for the O(v ) QCD amplitude are both UV and IR 
finite, while the 0(v 2 ) amplitude is UV finite albeit IR 
divergent. 

The occurrence of the IR divergences in the hard region 
at 0(a s v 2 ) is just as expected. From the pull-up mecha- 
nism, one can identify the IR divergences encountered in 
the hard-region calculation with those would arise from 
the soft (fc M ~ mv) region in a literal QCD-side calcu- 
lation, which must be canceled out upon matching. In 
fact, these IR divergences can be reconstructed with the 
knowledge of the 0(a s v 2 ) correction to the perturba- 
tive NRQCD matrix elements (cc( 3 5i)|V' t cr-ex|0)( 1 ) [23], 

(cc( 1 S* o )|V' t x|0) (1) [24], as well as c[, 0) [28]. Thus, we fi- 
nally end up with the both UV, IR-finite short-distance 
coefficients c| (i = 0,2). Our finding is consistent with 
the all-order-in-a s proof outlined in [38], that NRQCD 
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factorization holds for the exclusive production of one 
S'-wave quarkonium plus any higher orbital angular mo- 
mentum quarkonium in e + e~ annihilation. 

B. Asymptotic expressions of NLO short-distance 
coefficients 

The NRQCD short-distance coefficients Cg 1 ' and <£p 
are in general complex- valued. Their analytic expressions 

I 



■sH'4 



43 M. 4) 

' \ s m cJ 



(1) ( Mr M/ \ 
C 2.2 \ r ' i 2 



where (3q — ^-Ca ~ | n / is the one- loop coefficient of 
the QCD p function, and rij • = 4 denotes the number 
of active quark flavors. fi r denotes the renormalization 
scale, with the natural magnitude of order y/s; and [if is 
identified with the factorization scale in the MS scheme, 
whose value lies in somewhere between m c v and m c , the 
UV cutoff scale of NRQCD. 



It has been pointed out [39] that a peculiar double- 

2 (1) 

logarithmic correction oc In r arises in c for this 
helicity-flipped process, and our (9a) exactly agrees with 
the corresponding expression there. Equations (9) im- 
ply that the presence of the double logarithm persists to 
0{v 2 ) as well. Presumably, provided that y/s >• m c , one 
must resum these types of large logarithms to all orders 
to obtain a reliable prediction. Unfortunately, due to our 
incapability for the light-cone approach to dealing with 
the helicity-suppressed hard exclusive processes beyond 
LO in a s , currently it remains to be an open question 
how to fulfill such a resummation. 



are somewhat lengthy and will not be reproduced here. 
Nevertheless, it is enlightening to know their asymptotic 
expressions in the limit y/s » m c : 



I 

VI. PHENOMENOLOGY 

In the numerical analysis, we take \/s — 10.58 GeV, 
and the QED coupling constant a(\fs) = 1/130.9 [16]. 
The running strong coupling constant is evaluated by us- 
ing the two-loop formula with = 0.338 GeV [13, 14]. 
The LO NRQCD matrix elements are taken from [22]: 
(Oi)j/ip w (Oi)r, c = 0.387 GeV 3 . If we choose m c = 1.4 
GeV (r = 0.0700), the Gremm-Kapustin relation implies 
that (w 2 ).//^ = 0.223 and (v 2 ) Vc = 0.133. We will take 
[if = m c . 

We are ready to carry out a detailed analysis for the 
processes e + e~ — > J/ip + rj c and confront the B factory 
measurements. One important source of theoretical un- 
certainties comes from the scale setting for the strong 
coupling constant. There is no way to circumvent the 
scale ambiguity problem within the confine of NRQCD 
factorization, and we proceed to estimate the cross sec- 
tion by affiliating all the occurring a s with a common 
scale, [i r , and choosing [i r = y/s/2 and [i r = 2m c , respec- 
tively. It is hoped that the less biased results interpolate 
between these two sets of predictions. 

With r = 0.07, for [i r = y/s/2, we find c^/c^ = 



(o) \ n ( 1, s 5\ A3, 2 5, „, 41, 53, 2o 65 , „ 
= 4° ) x|/3„(- i ln^ + -j + (-lnV + i ln21n,--ln,--ln 2 2 + y ln2 

-^- 1 l) +i J 1 -P 0+ 1 llnr + 5 -ln2-^)\, (9 
36 4 J V 4 12 4 24/ J ' V 

1 (o) fl6, A n ( 1, s 11\ /3, 2 19, „, 31, 1, 2 
= ^ )x {y ln 4 +/30 ("4 ln 4^ + T2j + U ln V+ T2 ln21 ^ + 24 ln ^24 ln 

893, „ 5 2 497\ (\ n 3, 19 , „ 9\ 1 
+ ^ ta2 -^-^J +,,r U A,+ 4 lnr+ 12 In2+ 8j}' ( ° 

= 2 C (°) X [t In 4 +A)f-- In A + -) + (- 1^ r+- ln21nr-lln r - - In 2 2 
3 [3 ml MU V 4 4^ 2 3/ \12 12 24 8 

241 , „ 1 2 99 \ / 1 1 , 11 , „ 1 \1 
+ 144 ln2 -8 7r -T6j + ^U A+ 6 lnr+ T2 ln2 -24j[ (£ 
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TABLE I: Individual contributions to the predicted a[e + e~ — > 
J/i/> + r/c] at y/s = 10.58 GeV, labeled by powers of a s and v. 
The cross sections are in units of fb. 



a s (p r ) 


a (0) 


4 1 ' 


a (0) 

°2 


4 1] 


Q S (^) = 0.211 


4.40 


5.22 


1.72 


0.73 


a s (2m c ) = 0.267 


7.00 


7.34 


2.73 


0.24 



40 
30 
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3456789 10 
fir (GeV) 



FIG. 2: The /i r dependence of the cross section for e + e — > 
J/ip + r] c at y/s = 10.58 GeV. The 5 curves from bottom to top 
are erg ' (solid line), <Tq ' + o"2°' (dashed line), + OQ 1 ' (solid 
line), o-q^+ct^'+ctq 1 ' (dashed line), and 

(solid line), respectively. The blue and green bands represent 
the measured cross sections by the Belle and BaBar experi- 
ments, with respective systematic and statistical errors added 
in quadrature. 

8.83 - 6.02i, c£\/c£l = -5.09 + 8.89i, and c^/c^ = 
-5.15 + 8.99i: for pt r = 2m c we have /c[, 0) = 6.18 - 
6.02«, c^l/cf} = -7.74 + 8.89i, and c^/c^ = -7.80 + 

8.99i. Since Refcg 1 ]/^ -] is large and negative for i = 
1, 2, in conjunction with (7), one may think that the new 
0(a s v 2 ) correction would dilute the known 0{v 2 ) effect. 

Table I lists the predicted cross sections for e + e~ — > 
J/ip + rj c in double expansions of a s and u, with two 
sets of fi r . We reproduce the well-known results, i.e., the 
positive and significant 0(a s ) correction [13, 14], and the 
positive but less pronounced 0(v 2 ) correction [3, 15, 16]. 



The new 0(a s v 2 ) ingredient, cr^\ is positive but modest. 
It may be attributed to the near cancelation between the 
two terms in the prefactor of (v 2 )h in (6b). 

In Fig. 2, we plot the a[e + e~ — > J/ip + r] c ] as a func- 
tion of \i r . Whether including the contribution of or 
not clearly does not make a big difference. When [i r is 
relatively small, the state-of-the-art NRQCD prediction 
converges to the Babar measurement within errors. Had 
we taken somewhat larger values of the NRQCD matrix 
elements {Oi)h as in [15, 16], the agreement with the 
two B factories measurements would be better. 

VII. SUMMARY 

In this work we have computed the 0{a s v 2 ) correction 
to the helicity-suppressed process e + e~ — > J/ip + r\ c in 
the NRQCD factorization framework. The correspond- 
ing NLO perturbative short-distance coefficients associ- 
ated with the J/ip + ?y c EM form factor are directly ex- 
tracted from the hard loop-momentum region through 
relative order v 2 . By examining the asymptotic form of 
these coefficients, we reconfirm the pattern recognized in 
[22, 39]: The hard exclusive processes involving doublc- 
charmonium at higher twist in general are plagued with 
double logarithms of form In 2 s/m 2 once beyond LO in 
a s . When y/s 3> rn c , in order to obtain the reliable pre- 
dictions for such types of processes, one is enforced to 
resum these potentially large double logarithms to all or- 
ders in a s , which so far remains to be an open challenge. 

At the B factory energy, we found that incorporating 
this new piece of correction only modestly enhances the 
existing NRQCD predictions. This may be ascribed to 
some accidental cancelation between two different sources 
of relativistic corrections. At much higher y/s, this 
0{a s v 2 ) correction would be much more relevant. 
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